Assignment 2

Re-employment dataset, treatment 3 (№5)

This is the first time I am using an R project, and I apologize for submitting a plethora of files in the archive.

To install the written package, I use ‘build’ -> ‘more’ -> ‘load all.’ Again, I apologize for any clunkiness.

A word of introduction

The “exploring variability” section is commented out due to its lengthy runtime (takes an hour for each configuration). Similarly, the “combined treatment” part is also commented for brevity (but can be easily executed). The outcomes of these commented sections are detailed in the text tables.

Additionally, I have developed an R package called blueberry, containing all learners (except Casual Forest), along with visualizations and simulations. Brief documentation is provided for most of the functions.

1. Prepare and explore the data

I analyze the Pennsylvania re-employment bonus experiment (Bilias, 2000; Chernozhukov, 2018).

```r
Penn <- as.data.frame(read.table(\data/penn_jae.dat\, header = TRUE))

# To study all treatment arms simultaneously [--> COMMENT]: 
Penn    <- subset(Penn, tg==3 | tg==0)           # Keep only the 3rd treatment option
Penn$tg <- ifelse(Penn$tg == 3, 1, Penn$tg)      # make it {0,3} --> {0,1} 

# To study all treatment arms simultaneously [--> UNCOMMENT]: 
# Penn$tg <- ifelse(Penn$tg %in% c(2, 3, 4, 5, 6), 1, Penn$tg)

attach(Penn)

<!-- rnb-source-end -->

<!-- rnb-chunk-end -->


<!-- rnb-text-begin -->


<br>

These experiments were conducted in the 1980s by the U.S. Department of Labor to test the incentive effects of alternative compensation schemes for unemployment insurance (UI): UI claimants were randomly assigned either to a control group or one of six treatment groups.

I focus on **treatment group 3**, specifically, the one with a **high** bonus and **short** qualification:

| Group           | Bonus Amount             | Qualification Period | Workshop Offer |
|-----------------|--------------------------|----------------------|----------------|
| Control Group   | 0                        | 0                    | No             |
| Treatment 1     | Low (3 x WBA)            | Short (6 Weeks)      | Yes            |
| Treatment 2     | Low (3 x WBA)            | Long (12 Weeks)      | Yes            |
| **Treatment 3** | **High (6 x WBA)**       | **Short (6 Weeks)**  | **Yes**        |
| Treatment 4     | High (6 x WBA)           | Long (12 Weeks)      | Yes            |
| Treatment 5     | Initially High (6 x WBA) | Long (12 Weeks)      | Yes            |
| Treatment 6     | High (6 x WBA)           | Long (12 Weeks)      | No             |

: **Table 1**: Treatment Assignments

In the control group the current rules of the UI applied. Individuals in the treatment groups were offered a cash bonus if they found a job within some pre-specified period of time (qualification period), provided that the job was retained for a specified duration.

The matrix of covariates, $X$, consists of age group dummies, gender, race, the number of dependents, quarter of the experiment, location within the state, existence of recall expectations and type of occupation:

| Variable     | Description                                                                    |
|--------------|--------------------------------------------------------------------------------|
| `abdt`       | Chronological time of enrollment of each claimant in the experiment            |
| `tg`         | The treatment group (bonus amount - qualification period) of each claimant     |
| `inuidur1`   | A measure of the length (in weeks) of the first spell of unemployment          |
| `female`     | Dummy variable; it indicates if the claimant's sex is female (=1) or male (=0) |
| `black`      | Dummy variable; it indicates a person of black race (=1)                       |
| `hispanic`   | Dummy variable; it indicates a person of Hispanic race (=1)                    |
| `othrace`    | Dummy variable; it indicates a non-white, non-black, not-Hispanic (=1)         |
| `dep`        | The number of dependents of each claimant                                      |
| `q1-q6`      | Six dummy variables: the quarter of enrollment                                 |
| `recall`     | 1 if answered "yes" when asked if had any expectation to be recalled           |
| `agelt35`    | 1 if the claimant's age is less than 35                                        |
| `agegt54`    | 1 if the claimant's age is more than 54                                        |
| `durable`    | 1 if the occupation of the claimant was in the sector of durable manufacturing |
| `nondurable` | 1 if the claimant worked in the sector of non-durable manufacturing            |
| `lusd`       | 1 if filed in areas of low unemp. rate and short duration of unemployment      |
| `husd`       | 1 if filed in areas of high unemp. rate and short duration of unemployment     |
| `muld`       | 1 if filed in areas of moderate unemp. rate and long duration of unemployment  |

<br>

## 2. Descriptive statistics & simple ATE

The data comes from an experimental (RCT) setting:


<!-- rnb-text-end -->


<!-- rnb-chunk-begin -->


<!-- rnb-source-begin eyJkYXRhIjoiYGBgclxuYGBgclxuaGVhZChQZW5uKVxuYGBgXG5gYGAifQ== -->

```r
```r
head(Penn)

<!-- rnb-source-end -->

<!-- rnb-frame-begin eyJtZXRhZGF0YSI6eyJjbGFzc2VzIjoiZGF0YS5mcmFtZSIsIm5yb3ciOjYsIm5jb2wiOjIzLCJzdW1tYXJ5Ijp7IkRlc2NyaXB0aW9uIjoiZGYgWzYgw5cgMjNdIn19LCJyZGYiOiJINHNJQUFBQUFBQUFBN1ZVelU3RE1BeDIybTVqUlV3VE8zQ0FGd0Jway9hSHVIRkJpRE1DYWRlczdicHFXVnJhRkhncW5vZkg0STQwY0x1a3k2cU5IdzBpV2Y3c3o0NGRwK250MWFodmoyd0FNTUV5RFRBckNLRnlmM2ZkdmdDd0REUUlXRkRQOURNR3RSQWNvUnlnVkFIT2JwWnkrbzd5Z3JpSitoVzV4cExmdkM3Zk51dlZ2bkNZRndZNHpqc0RPSkgrYmZ4K2lkY1hrUUpiZUgzOUZXOW9Qa1BULzFHZmZNRVRUUk9OSjZXOFhlcnYwdDl2OWlmZjhPWDhUV2ZVei8vVCttUWxhdytpd3VuY1N4QTBJWHNVUzZkRng2NlEyQkMrUkhzQlR3TTNqYnNsdXlmdDZzU2JVK2FwamNlTU9qTVZPZzJTaVBMQWtYWXRGTk9ZT2lyV2RMMUlsWHZvRnFoWG9INkJCZ1VhRnVoY05SQjdEbVZNMWFDK3gwUi9xSm0rR0tyOEdqWk94MFczTmcvNXVzZGlhZUlxUE5Yd1BHVnVhWWoxT0h6cXFFSHF0Mk5sZzhqT253VkIvc0xYcCs4d21xanBLNmZ0VWtFN2t4ajNRMnRSU3FtRmtRaENqa2xHUzFiU2swbGNjalJUbm5YbXRwMXB5bWZ0N3FERW14SDNzNkt3L3I0SXJQNUZKUDgwOGphTUQ1bGVWVFAzdUIvdzR0SVpIWHRzZGF1UEVqWndRUGw4T2xFY2NQVnAyZWhOT2lJVVZLWFlUc2lVSno4NkxENEJabFdRN3QwRkFBQT0ifQ== -->

<div data-pagedtable="false">
  <script data-pagedtable-source type="application/json">
{"columns":[{"label":[""],"name":["_rn_"],"type":[""],"align":["left"]},{"label":["abdt"],"name":[1],"type":["int"],"align":["right"]},{"label":["tg"],"name":[2],"type":["dbl"],"align":["right"]},{"label":["inuidur1"],"name":[3],"type":["int"],"align":["right"]},{"label":["inuidur2"],"name":[4],"type":["int"],"align":["right"]},{"label":["female"],"name":[5],"type":["int"],"align":["right"]},{"label":["black"],"name":[6],"type":["int"],"align":["right"]},{"label":["hispanic"],"name":[7],"type":["int"],"align":["right"]},{"label":["othrace"],"name":[8],"type":["int"],"align":["right"]},{"label":["dep"],"name":[9],"type":["int"],"align":["right"]},{"label":["q1"],"name":[10],"type":["int"],"align":["right"]},{"label":["q2"],"name":[11],"type":["int"],"align":["right"]},{"label":["q3"],"name":[12],"type":["int"],"align":["right"]},{"label":["q4"],"name":[13],"type":["int"],"align":["right"]},{"label":["q5"],"name":[14],"type":["int"],"align":["right"]},{"label":["q6"],"name":[15],"type":["int"],"align":["right"]},{"label":["recall"],"name":[16],"type":["int"],"align":["right"]},{"label":["agelt35"],"name":[17],"type":["int"],"align":["right"]},{"label":["agegt54"],"name":[18],"type":["int"],"align":["right"]},{"label":["durable"],"name":[19],"type":["int"],"align":["right"]},{"label":["nondurable"],"name":[20],"type":["int"],"align":["right"]},{"label":["lusd"],"name":[21],"type":["int"],"align":["right"]},{"label":["husd"],"name":[22],"type":["int"],"align":["right"]},{"label":["muld"],"name":[23],"type":["int"],"align":["right"]}],"data":[{"1":"10824","2":"0","3":"18","4":"18","5":"0","6":"0","7":"0","8":"0","9":"2","10":"0","11":"0","12":"0","13":"0","14":"1","15":"0","16":"0","17":"0","18":"0","19":"0","20":"0","21":"0","22":"1","23":"0","_rn_":"1"},{"1":"10824","2":"0","3":"1","4":"1","5":"0","6":"0","7":"0","8":"0","9":"0","10":"0","11":"0","12":"0","13":"0","14":"1","15":"0","16":"0","17":"0","18":"0","19":"0","20":"0","21":"1","22":"0","23":"0","_rn_":"4"},{"1":"10747","2":"0","3":"27","4":"27","5":"0","6":"0","7":"0","8":"0","9":"0","10":"0","11":"0","12":"0","13":"1","14":"0","15":"0","16":"0","17":"0","18":"0","19":"0","20":"0","21":"1","22":"0","23":"0","_rn_":"5"},{"1":"10670","2":"1","3":"3","4":"3","5":"1","6":"0","7":"0","8":"0","9":"2","10":"0","11":"0","12":"1","13":"0","14":"0","15":"0","16":"0","17":"1","18":"0","19":"0","20":"0","21":"0","22":"0","23":"1","_rn_":"8"},{"1":"10768","2":"1","3":"28","4":"11","5":"1","6":"0","7":"0","8":"0","9":"0","10":"0","11":"0","12":"0","13":"1","14":"0","15":"0","16":"0","17":"0","18":"0","19":"0","20":"0","21":"0","22":"0","23":"1","_rn_":"9"},{"1":"10712","2":"1","3":"6","4":"6","5":"0","6":"0","7":"0","8":"0","9":"2","10":"0","11":"0","12":"0","13":"1","14":"0","15":"0","16":"0","17":"0","18":"0","19":"1","20":"0","21":"0","22":"0","23":"1","_rn_":"11"}],"options":{"columns":{"min":{},"max":[10],"total":[23]},"rows":{"min":[10],"max":[10],"total":[6]},"pages":{}}}
  </script>
</div>

<!-- rnb-frame-end -->

<!-- rnb-chunk-end -->


<!-- rnb-text-begin -->


Interestingly, there are no continuous variables apart from the response variable.

<br>

### Covariate balance check

This is done with the Eicher-Huber-White (heteroscedasticity adjusted) standard errors, including the linear terms as well as the second order interactions:


<!-- rnb-text-end -->


<!-- rnb-chunk-begin -->


<!-- rnb-source-begin eyJkYXRhIjoiYGBgclxuYGBgclxubSA8LSBsbSh0Z34oZmVtYWxlK2JsYWNrK290aHJhY2UrZmFjdG9yKGRlcCkrcTIrcTMrcTQrcTUrcTYrYWdlbHQzNSthZ2VndDU0K2R1cmFibGUrbHVzZCtodXNkKV4yKVxubSA8LSBjb2VmdGVzdChtLCB2Y292ID0gdmNvdkhDKG0sIHR5cGU9XFxIQzFcXCkpXG5cbnN0YXJnYXplcihtWzE6MTUsXSwgdHlwZSA9IFxcdGV4dFxcLHRpdGxlID0gJ0NvdmFyaWF0ZSBiYWxhbmNlIGNoZWNrJylcbmBgYFxuYGBgIn0= -->

```r
```r
m <- lm(tg~(female+black+othrace+factor(dep)+q2+q3+q4+q5+q6+agelt35+agegt54+durable+lusd+husd)^2)
m <- coeftest(m, vcov = vcovHC(m, type=\HC1\))

stargazer(m[1:15,], type = \text\,title = 'Covariate balance check')

<!-- rnb-source-end -->

<!-- rnb-output-begin eyJkYXRhIjoiXG5Db3ZhcmlhdGUgYmFsYW5jZSBjaGVja1xuPT09PT09PT09PT09PT09PT09PT09PT09PT09PT09PT09PT09PT09PT09PT09PT09PT09PVxuICAgICAgICAgICAgIEVzdGltYXRlIFN0ZC4gRXJyb3IgdCB2YWx1ZSBQcig+IHwgdHwgKVxuLS0tLS0tLS0tLS0tLS0tLS0tLS0tLS0tLS0tLS0tLS0tLS0tLS0tLS0tLS0tLS0tLS0tLVxuKEludGVyY2VwdCkgICAwLjIwOCAgICAgMC4xNTUgICAgIDEuMzQ1ICAgICAwLjE3OSAgIFxuZmVtYWxlICAgICAgICAtMC4xMTkgICAgMC4xMzYgICAgLTAuODc2ICAgICAwLjM4MSAgIFxuYmxhY2sgICAgICAgICAwLjE2MyAgICAgMC4wNzkgICAgIDIuMDU4ICAgICAwLjA0MCAgIFxub3RocmFjZSAgICAgICAwLjA0NCAgICAgMC4yODcgICAgIDAuMTUyICAgICAwLjg3OSAgIFxuZmFjdG9yKGRlcCkxICAwLjQzMCAgICAgMC4xNzUgICAgIDIuNDU0ICAgICAwLjAxNCAgIFxuZmFjdG9yKGRlcCkyICAwLjEzMCAgICAgMC4xNTIgICAgIDAuODUzICAgICAwLjM5NCAgIFxucTIgICAgICAgICAgICAwLjE0MSAgICAgMC4xNTYgICAgIDAuOTAzICAgICAwLjM2NyAgIFxucTMgICAgICAgICAgICAwLjExNSAgICAgMC4xNTYgICAgIDAuNzM5ICAgICAwLjQ2MCAgIFxucTQgICAgICAgICAgICAwLjEzMyAgICAgMC4xNTYgICAgIDAuODU2ICAgICAwLjM5MiAgIFxucTUgICAgICAgICAgICAwLjEwNiAgICAgMC4xNTUgICAgIDAuNjgzICAgICAwLjQ5NCAgIFxucTYgICAgICAgICAgICAwLjAzMiAgICAgMC4xNTQgICAgIDAuMjA5ICAgICAwLjgzNCAgIFxuYWdlbHQzNSAgICAgICAwLjEwMyAgICAgMC4xMzAgICAgIDAuNzkzICAgICAwLjQyOCAgIFxuYWdlZ3Q1NCAgICAgICAwLjM0NiAgICAgMC4yMjMgICAgIDEuNTUyICAgICAwLjEyMSAgIFxuZHVyYWJsZSAgICAgICAtMC4xMzIgICAgMC4xNjcgICAgLTAuNzkwICAgICAwLjQzMCAgIFxubHVzZCAgICAgICAgICAwLjE0MCAgICAgMC4wNjggICAgIDIuMDY1ICAgICAwLjAzOSAgIFxuLS0tLS0tLS0tLS0tLS0tLS0tLS0tLS0tLS0tLS0tLS0tLS0tLS0tLS0tLS0tLS0tLS0tLVxuIn0= -->

Covariate balance check

Estimate Std. Error t value Pr(> | t| )

(Intercept) 0.208 0.155 1.345 0.179
female -0.119 0.136 -0.876 0.381
black 0.163 0.079 2.058 0.040
othrace 0.044 0.287 0.152 0.879
factor(dep)1 0.430 0.175 2.454 0.014
factor(dep)2 0.130 0.152 0.853 0.394
q2 0.141 0.156 0.903 0.367
q3 0.115 0.156 0.739 0.460
q4 0.133 0.156 0.856 0.392
q5 0.106 0.155 0.683 0.494
q6 0.032 0.154 0.209 0.834
agelt35 0.103 0.130 0.793 0.428
agegt54 0.346 0.223 1.552 0.121
durable -0.132 0.167 -0.790 0.430
lusd 0.140 0.068 2.065 0.039
—————————————————-

<!-- rnb-output-end -->

<!-- rnb-chunk-end -->


<!-- rnb-text-begin -->


Only a portion is printed for brevity. Generally, the covariates are balanced, with very few exceptions. According to the F-test, the covariates are not jointly significant at all conventional significance levels:

F-statistic: 0.9825 on 102 and 5136 DF, p-value: 0.5315

I assume that the RCT assumptions hold.

<br>

### Simple ATE estimation

To evaluate the impact of the treatment on unemployment duration, I consider the linear regression model:

$$ Y = W\beta_1 + \mathbf{X}'\beta_2 + \epsilon, \quad E(\epsilon | W, \mathbf{X}')' = 0, $$

where $Y$ is the **log** of the duration of unemployment, $W$ is a treatment indicator, and $\mathbf{X}$ is a set of controls defined previously. Here $\beta_1$ is the ATE, if the RCT assumptions hold.

I also consider an interactive regression model:

$$ Y = W\alpha_1 + W \mathbf{X}'\alpha_2 + \mathbf{X}'\beta_2 + \epsilon, \quad E(\epsilon | W, \mathbf{X}', W\mathbf{X}')' = 0, $$

where $\mathbf{X}'$'s are demeaned (apart from the intercept), so that $\alpha_1$ is the ATE, if the RCT assumptions hold.

Under RCT, the projection coefficient $\beta_1$ has the interpretation of the causal effect of the treatment on the average outcome (ATE). The covariates here are independent of the treatment $W$, so we can identify $\beta_1$ by just linear regression of $Y$ on $W$, without adding covariates. However, I do add covariates to improve the precision of the ATE estimates.

I specify four models:

1.  Simple treated-control mean comparison,

2.  Adding covariates as controls,

3.  Running an interactive regression model (*including all variables and second-order interactions*),

4.  Partialling out with Lasso.

Last two methods are taken from Viktor Chernozukov's [kaggle](https://www.kaggle.com/code/victorchernozhukov/analyzing-rct-reemployment-experiment/notebook).

The omitted dummies, ensuring no perfect multicollinearity, are: `q1`, `nondurable`, `muld`.


<!-- rnb-text-end -->


<!-- rnb-chunk-begin -->


<!-- rnb-source-begin eyJkYXRhIjoiYGBgclxuIyAxXSBOYWl2ZSBtZWFuIGNvbXBhcmlzb24gXG5vbHMuc2ltcGxlICA9ICBsbShsb2coaW51aWR1cjEpfnRnKVxub2xzLnNpbXBsZSAgPSAgY29lZnRlc3Qob2xzLnNpbXBsZSwgdmNvdiA9IHZjb3ZIQyhvbHMuc2ltcGxlLCB0eXBlPVwiSEMxXCIpKVxuXG4jIDJdIFJlZ3Jlc3Npb24gd2l0aCBjb250cm9sc1xub2xzLmNvbnRyb2wgPSBsbShsb2coaW51aWR1cjEpfnRnKyAoZmVtYWxlK2JsYWNrK290aHJhY2UrZmFjdG9yKGRlcCkrcTIrcTMrcTQrcTUrcTYrYWdlbHQzNSthZ2VndDU0K2R1cmFibGUrbHVzZCtodXNkKV4yKVxub2xzLmNvbnRyb2wgPSBjb2VmdGVzdChvbHMuY29udHJvbCwgdmNvdiA9IHZjb3ZIQyhvbHMuY29udHJvbCwgdHlwZT1cIkhDMVwiKSlcblxueCA9IG1vZGVsLm1hdHJpeCAofiAoZmVtYWxlK2JsYWNrK290aHJhY2UrZmFjdG9yKGRlcCkrcTIrcTMrcTQrcTUrcTYrYWdlbHQzNSthZ2VndDU0K2R1cmFibGUrbHVzZCtodXNkKV4yKVssLTFdXG54ID0gc2NhbGUoeCwgY2VudGVyID0gVFJVRSwgc2NhbGUgPSBGQUxTRSlcblxuIyAzXSBQYXJ0aWFsbGluZyBvdXQgd2l0aCBsYXNzb1xudGdfZGVtZWFuICA9IHRnIC0gbWVhbih0ZylcbkR4ICAgICAgICAgPSBtb2RlbC5tYXRyaXgofnRnX2RlbWVhbip4KVssLTFdXG5ybGFzc28uaXJhID0gc3VtbWFyeShybGFzc29FZmZlY3RzKER4LCBsb2coaW51aWR1cjEpLCBpbmRleCA9IDEpKVxuXG4jIDRdIEludGVyYWN0aXZlIHJlZ3Jlc3Npb24gbW9kZWwgXG5vbHMuaXJhID0gbG0obG9nKGludWlkdXIxKSB+IHRnKngpIFxub2xzLmlyYSA9IGNvZWZ0ZXN0KG9scy5pcmEsIHZjb3YgPSB2Y292SEMob2xzLmlyYSwgdHlwZT1cIkhDMVwiKSlcblxub2xzLmlyYV9vcmQ8LSBvbHMuaXJhW29yZGVyKG9scy5pcmFbLCBcIlByKD58dHwpXCJdKSwgXVxuc3RhcmdhemVyKG9scy5pcmFfb3JkWzE6MTUsXSwgdHlwZSA9IFwidGV4dFwiLHRpdGxlID0gJ0ludGVyYWN0aXZlIHJlZ3Jlc3Npb24gbW9kZWwnKVxuYGBgIn0= -->

```r
# 1] Naive mean comparison 
ols.simple  =  lm(log(inuidur1)~tg)
ols.simple  =  coeftest(ols.simple, vcov = vcovHC(ols.simple, type="HC1"))

# 2] Regression with controls
ols.control = lm(log(inuidur1)~tg+ (female+black+othrace+factor(dep)+q2+q3+q4+q5+q6+agelt35+agegt54+durable+lusd+husd)^2)
ols.control = coeftest(ols.control, vcov = vcovHC(ols.control, type="HC1"))

x = model.matrix (~ (female+black+othrace+factor(dep)+q2+q3+q4+q5+q6+agelt35+agegt54+durable+lusd+husd)^2)[,-1]
x = scale(x, center = TRUE, scale = FALSE)

# 3] Partialling out with lasso
tg_demean  = tg - mean(tg)
Dx         = model.matrix(~tg_demean*x)[,-1]
rlasso.ira = summary(rlassoEffects(Dx, log(inuidur1), index = 1))

# 4] Interactive regression model 
ols.ira = lm(log(inuidur1) ~ tg*x) 
ols.ira = coeftest(ols.ira, vcov = vcovHC(ols.ira, type="HC1"))

ols.ira_ord<- ols.ira[order(ols.ira[, "Pr(>|t|)"]), ]
stargazer(ols.ira_ord[1:15,], type = "text",title = 'Interactive regression model')

Sorted by p-value. Only a part is printed for brevity. Based on the IRA results one may conclude that there is a possible treatment effect heterogeneity: some interactions between \(X\) and \(W\) are significant.

OLS_results(ols.simple, ols.control, ols.ira, rlasso.ira)

Treatment group №3 experiences an average decrease of about 5.8% in the length of unemployment spell, albeit the effects are not significant at the 5% level.

We also see the regression estimators offer slightly lower estimates: the difference is perhaps due to minor imbalance in the treatment allocation, which the regression estimators try to correct.

Study’s report conclusion:
Estimates based on data from UI wage records fail to provide any evidence that the bonus offers enhanced the employment of claimants following the benefit application date.

3. Modeling choices and heterogeneity detection

The choice of propensity score

Firstly, I define the standard variables necessary for future estimation (\(Y,X,W\)):

# Produce matrix of covariates X, treatment indicator W and outcome vector Y: 
X <- model.matrix(~ (female+black+othrace+factor(dep)+q2+q3+q4+q5+q6+agelt35+agegt54+durable+lusd+husd))[,-1]
W <- tg
Y <- log(inuidur1)

In the third treatment arm, there are 1,885 treated units and 3,354 units in the control group:

table(tg)

As the treatment is randomly assigned, I use the fraction of treated as the estimator of the propensity score (follows Chernozhukov et al.,2018). Hence, the propensity score is constant: this value is used by learners for the estimation of CATE.

# e_rct = rep(e_rct, length(Y))
e_rct = mean(W)
e_rct

To double check: propensity scores estimated by regression forest are tightly clustered around mean(W).

# Propensity model
W_hat <- regression_forest(X, W, tune.parameters = "all")
W_hat.rf  <- W_hat$predictions

hist(W_hat.rf, main = "Estimated propensity scores", xlab = "Propensity Score", xlim = c(0, 1))

Heterogeneity detection

I begin by testing pre-specified null hypotheses of the form

\[ H_{0}: \mathop{\mathrm{E}}[Y(1) - Y(0) | G_i = 1] = \mathop{\mathrm{E}}[Y(1) - Y(0) | G_i = 0]. \]

That is, that the treatment effect is the same regardless of membership to the group \(G_i\). I choose to examine how the outcome of participating in the experiment differs for younger people (agelt35 = 1).

When predicting the conditional average treatment effect \(\tau(X_i)\) for observation \(i\) we should in general avoid using a model that was fitted using observation \(i\) This sort of sample splitting (honesty) is one of the required ingredients to get unbiased estimates of the CATE. I do 3 fold cross-fitting:

Y_tilde = aipw_GATE(X,W,Y,e_rct,3)
summary(lm_robust(Y_tilde ~ X[, 'agelt35']))

We observe that younger people may experience approximately a 14% reduction in the duration of unemployment spells compared to the average program effect. However, this result is not robust across different seeds and numbers of folds.

With no further hypotheses to explore, I provide a summary of the underlying effect heterogeneity:

summary(lm_robust(Y_tilde ~ X))

The coefficients in this regression suggest general trends, but they should not be interpreted as partial effects. Such interpretation would only be accurate if the true model were genuinely linear in covariates — an assumption we generally shouldn’t be willing to make.

While some coefficients in the interactive regression may appear significant, they could be false positives. A more aggregated AIPW regression indicates no linear effect heterogeneity.

4. Split the data into training and test sample

set.seed(123)  

indices   <- sample(1:nrow(X), size = 0.5*nrow(X))
train_idx <- numeric(nrow(Penn))

train_idx[indices]  <- 1  

X_train = X[indices,]
X_test  = X[-indices,]

Y_train = Y[indices]
Y_test  = Y[-indices]

W_train = W[indices]
W_test  = W[-indices]

# Summary of treated and control in train and test samples:
table(Train = ifelse(Penn[train_idx==1, ]$tg == 0, "Control", "Treated"))
table(Test = ifelse(Penn[train_idx==0, ]$tg == 0, "Control", "Treated"))

5. Apply Causal ML methods

Before applying any machine learning methods, I think it is important to discuss the estimation methods. The following reasoning is based on Okasa (2022).

Their simulation evidence suggests that using the full sample for estimation of both the nuisance functions as well as the CATE function leads to a remarkably good performance in terms of both bias and variance in finite samples.

The possible reason for this phenomenon might be due to the out-of-bag predictions of the forest that are used throughout the simulations. Even though these predictions are not out-of-sample per se they are not directly based on the observations used for estimation and as such might help to alleviate the overfitting problem when using full sample.

Using cross-fitted estimation, one should (theoretically) observe a smaller bias but higher variance of the estimators. However, in almost all cases considered in the paper the authors observe both higher bias as well as higher variance, particularly for the small sample sizes.

Based on the above simulation evidence, it seems reasonable to always use the full-sample estimation together with out-of-bag predictions when a relatively small sample is available, whereas to use the cross-fitting procedure when a relatively large data is accessible, regardless of the choice of a meta-learner.

Given that the total number of observations in the Penn dataset is limited (\(n=5239\)), I prefer to do the full-sample estimation.

S-learner

1. Use ML estimator to fit model \(\hat{m}(W,X)\) in the whole training sample

2. Estimate CATE as \(\hat{\tau}(x) = \hat{m}(1,x) - \hat{m}(0,x)\)

If the treatment indicator is not strongly predictive for the outcome, the S-learner will tend to estimate a zero treatment effect.

To verify \(W\)’s low predictive power, I first estimate whether post-lasso chooses the treatment variable \(W\) when predicting \(Y\):

post_lasso = hdm::rlasso(x = cbind(W,X),y = Y)
post_lasso$coefficients[post_lasso$coefficients != 0]

The treatment \(W\) is not selected. I then do the similar operation but using a decision tree:

fit <- rpart(Y~ W+X,control = rpart.control(cp = 0.0015))
rpart.plot(fit)

For this dataset, I conclude that the treatment \(W\) is not a strong predictor for the outcome \(Y\), which could potentially undermine the performance of the S-learner. Moreover, based on the contemporary literature:

1. Okasa (2022) does not recommend the usage of the S-learner for estimation of heterogeneous causal effects due to empirically documented undesirable statistical properties such as high bias and consistency issues.
2. Alaa and van der Schaar (2018) similarly find the T-learner to uniformly outperform the S-learner in estimating heterogeneous treatment effects.
3. Related simulation studies (Jacob, 2020; Knaus et al., 2021) find also a relatively competitive performance of the T-learner, especially with the Random Forest as a base learner.

Hence, I have decided to exclude this method due to unfavorable estimation conditions.

T-Learner

For the T-learner, I implement the following procedure:

1. Use regression forest to fit model \(\hat{m}(1,X)\) in treated

2. Use regression forest to fit model \(\hat{m}(0,X)\) in controls

3. Estimate CATE as \(\hat{\tau}(x) = \hat{m}(1,x) - \hat{m}(0,x)\)

This procedure is expected to work particularly well if the CATE function is complicated and there are no common trends in the response functions. It is expected to work rather poorly if the CATE function is simple, as the response functions are not trained jointly and their difference is unstable (Okasa, 2022).

No additional sample-splitting induced by multiple nuisance functions is required as the response functions are themselves estimated on separate samples defined by treated and control.

cate_tl = T_learner(X,W,Y,train_idx)
hist(cate_tl)

X-learner

The performance of the T-learner is aggravated if the treatment assignment is highly unbalanced, i.e. when the majority of observations in the sample belong to only one treatment status.

Since there are multiple treatment arms, the numbers of treated and control units are unbalanced, with 1,885 treated units and 3,354 units in control. Hence, we should expect some form of regularization bias for \(\hat{m}(1,X)\).

The following implementation is based on a Short Course by Athey, Spiess and Wager.

The first estimation step is identical to the T-learner. I use regression forest to fit model \(\hat{m}(1,X)\) in treated and \(\hat{m}(0,X)\) in control: \[\begin{align*} \hat{m}_0(X) = E[Y| W=0, X] \\ \hat{m}_1(X) = E[Y| W=1, X] \\ \end{align*}\]

For the second stage, I impute the treatment effect for the control and for the treated using: \[ \begin{align*} \hat{\tau}(x, 0) &= \hat{m}_1(x, 0) - Y \\ \hat{\tau}(x, 1) &= Y - \hat{m}_0(x, 1) \\\end{align*} \]

Lastly, I fit two more models to predict these effects. With \(\hat{e}(x)\) being the propensity score model, one can combine the two second stage models as follows:\[\begin{align*} \hat{\tau}(x) &= \hat{\tau}(x, 0) \hat{e}(x) + \hat{\tau}(x, 1) (1-\hat{e}(x)) \end{align*}\]

X-learner is expected to work particularly well in unbalanced settings, which is often the case in practice as the share of treated might be restricted financially or otherwise (Okasa, 2022).

Moreover, the X-learner might be less prone to overfitting bias which would partly justify the full-sample estimation as described in Künzel et al. (2019).

cate_xl = X_learner(X,W,Y,train_idx,e_rct)
hist(cate_xl)

Causal Forest

The causal forest estimates CATEs using a localized version of the partially linear estimator:

\[\begin{equation} \hat{\tau}^{cf}(x) = argmin_{\tau} \left\{ \sum_{i=1}^N \alpha_i(x) \left[(Y_i - \hat{m}(X_i)) - \tau(x)(W_i - \hat{e}(X_i)) \right]^2 \right\}, \end{equation} \] where \(\alpha(x)\) are \(x\)-specific weights.

Using predict function, I get out-of-bag estimated individual ATEs for every individual in the sample:

cf = causal_forest(X_train, Y_train, W_train, W.hat = e_rct)
cate_cf = predict(cf,X_test)$predictions
hist(cate_cf)

The grf package also produces a measure of variable importance that indicates how often a variable was used in a tree split.

This can be a rough diagnostic, but it should not be interpreted as indicating that, for example, variable with low importance is not related to heterogeneity: if two covariates are highly correlated, the trees might split on one covariate but not the other, even though both (or maybe neither) are relevant in the true DGP.

covariates <- c("female", "black", "othrace", "factor(dep)1", "factor(dep)2",
                     "q2", "q3", "q4", "q5", "q6", "agelt35", "agegt54", "durable",
                     "lusd", "husd")

var_imp        <- c(variable_importance(cf))
names(var_imp) <- covariates
sorted_var_imp <- sort(var_imp, decreasing = TRUE)

for (i in 1:5) {cat(sprintf("%-15s %-10.4f\n", names(sorted_var_imp)[i], sorted_var_imp[i]))}

DR-learner

The DR-learner creates the pseudo-outcome in a first step \[ \begin{align} \tilde{Y}_{ATE} = \underbrace{\hat{m}(1,X) - \hat{m}(0,X)}_{\text{outcome predictions}} + \underbrace{\frac{W (Y - \hat{m}(1,X))}{\hat{e}(X)} - \frac{(1-W) (Y - \hat{m}(0,X))}{1-\hat{e}(X)}}_{\text{weighted residuals}} \nonumber \end{align}\]and uses it in a final regression step (a generic ML problem).

DR-learner is expected to work well in similar situations as the X-learner with a more balanced treatment assignment, as too extreme propensity scores might yield the estimator unstable (not an issue for Penn). It should have an additional advantage over the X-learner due to its double robustness property.

• The simulations of Kennedy (2020) suggest a faster convergence rate of the DR-learner in comparison to the X- and T-learner.

• For the smaller sample sizes (500 and 2′000), the reduction in the overfitting bias is not large enough in comparison to the bias stemming from the estimation of the CATE function. For small sample sizes the additional splitting of the sample does not leave enough observations to learn the non-linear structure of the CATE (Okasa 2022).

It is well-suited for the nuisance parameters, but should have a harder time with the linear heterogeneous effects:

cate_dr = DR_learner(X,W,Y,train_idx,e_rct)
hist(cate_dr)

R-learner

The generic version of R-learner is to use the unmodified covariates in a weighted regression with pseudo outcomes: \[\hat{\tau}^{rl}(x) = argmin_{\tau} \sum_{i=1}^N \underbrace{(W_i - \hat{e}(X_i))^2}_{\text{weight}} \left(\underbrace{\frac{Y_i - \hat{m}(X_i)}{W_i - \hat{e}(X_i)}}_{\text{pseudo-outcome}} - X_i\beta \right)^2\]

R-learner is expected to work in a settings where the nuisance functions as well as the CATE function might have a high degree of complexity. Nie and Wager (2021) show good performance of the R-learner in settings with complicated nuisance functions and rather simple CATE function.

The performance of the R-learner is competitive with the other meta-learners especially in smaller samples, particularly for the full-sample version (Okasa, 2022).

cate_rl_rf = R_learner(X,W,Y,train_idx,e_rct)
hist(cate_rl_rf)

6. Predict the CATEs in the test sample

Now, let’s compare the estimated effects of the different methods used so far:

df <- data.frame(
  T_learner = cate_tl,
  X_learner = cate_xl,
  R_learner = cate_rl_rf,
  DR_learner = cate_dr,
  Causal_Forest = cate_cf)

df_long <- tidyr::gather(df, key = "Method", value = "CATE")

ggplot(df_long, aes(x = CATE, fill = Method)) +
  geom_density(alpha = 0.5) +
  labs(title = "Densities of Estimated CATE",
       x = "Conditional Average Treatment Effect (CATE)",
       y = "Density") + theme_minimal() + xlim(c(-1.5, 1.5)) # + ylim(c(0,3.2))

results = cbind(cate_tl,cate_xl,cate_cf,cate_dr,cate_rl_rf)
colnames(results) = c("T-learner","X-learner","Causal Forest","DR-learner","R-learner")
pairs.panels(results,method = "pearson")
describe(results)

All methods are highly correlated. DR and R-learners are more spread out in the predicted CATE. At first glance, the results look promising, without unexpected glitches.

7. Evaluating Heterogeneous Effects

One might be tempted to think: “if the histogram is concentrated at a point, then there is no heterogeneity; if the histogram is spread out, then our estimator has found interesting heterogeneity.” However, this may be false.

If the histogram is concentrated at a point, we may be under-powered: our method was not able to detect any heterogeneity, but maybe it would detect it if we had more data. If the histogram is spread out, we may be overfitting: our model is producing very noisy estimates \(\hat{\tau}(x)\), but in fact the true \({\tau}(x)\) can be much smoother as a function of \(x\).

The goal is to try to assess the validity of the estimated CATEs.

By hand: Best Linear Predictor

Strategy A: Weighted residual BLP

\[(\beta_1, \beta_2) = argmin_{b_1,b_2}E\{\omega(X)[Y - b_1 (W - e(X)) - b_2 (W - e(X)) (\hat{\tau}(X) - E[\hat{\tau}(X)]) - a \tilde{X}]^2\} \]

where \(\omega(X) = [e(X) (1-e(X))]^{-1}\).

Strategy B Horvitz-Thompson BLP\[(\beta_1, \beta_2) = argmin_{b_1,b_2}E\{[HY - b_1 - b_2 (\hat{\tau}(X) - E[\hat{\tau}(X)]) - a H\tilde{X}]^2\} \]

where \(H = \frac{W-e(X)}{e(X)(1-e(X))}\). The de-noising component \(\tilde{X}\) in both strategies is \(\hat{m}(0,X)\) obtained by regression forest.

\(\beta_2 = 0\) if either the effects are homogeneous or our CATE estimates are “bad.” Thus, rejecting the hypothesis \(\beta_2 = 0\) would imply that both the effects are heterogeneous and our CATE estimates are reliable - evaluCATE documentation

To compare all methods, I begin by using a modified function introduced in the Exercise 8. This function takes predicted CATE, test outcomes, test treatment, and de-noising variables as inputs, producing the estimated \(\beta_2\) coefficients and their standard errors as outputs.

# Denoising 
rfm0     = regression_forest(X_train[W_train==0,],Y_train[W_train==0])
mhat0_rf = predict(rfm0, X_test)$predictions

BLP_tl = BLP_lecture(cate_tl,Y_test,W_test,mhat0_rf,"T-learner RF")
BLP_xl = BLP_lecture(cate_xl,Y_test,W_test,mhat0_rf,"X-learner RF")
BLP_cf = BLP_lecture(cate_cf,Y_test,W_test,mhat0_rf,"Causal Forest")
BLP_dr = BLP_lecture(cate_dr,Y_test,W_test,mhat0_rf,"DR-learner")
BLP_rl = BLP_lecture(cate_rl_rf,Y_test,W_test,mhat0_rf,"R-learner")

results_BLP = list(BLP_tl,BLP_xl,BLP_cf,BLP_dr,BLP_rl)
names       = c("T-learner","X-learner","Causal Forest","DR-learner","R-learner")

plot_BLP(results_BLP, names)

It is worth mentioning that de-noising helps to get more precise estimates. Nonetheless, the results are quite miserable. Causal Forest predictions seem to be the best candidate for the evaluCATE.

Using evaluCATE package

“It’s like putting the real thing on the left hand side, but we don’t observe the real thing, but can approximate it using an unbiased signal \(HY\). If our predictions correlate with this unbiased signal, then we did a good job”, - M.C. Knaus

# Full sample (train+test) prediction: 
eval_cate_cf = predict(cf,X)$predictions

# Logical treatment indicator: 
train_idx_l  = as.logical(train_idx)

evaluation <- evaluCATE(Y, W, X, eval_cate_cf, train_idx_l, pscore = rep(e_rct, length(Y)))
summary(evaluation, target = "BLP")

The results are as follows:

• ATE \(\approx - 0.05\), but never statistically significant.

• HTE \(\approx 0.7\), and very noisy. We are partially able to detect the heterogeneity (and probably slightly overfit the CATE function: \(\beta_2 < 1\)).

The slope \(\beta_2\) is a measure of how the CATE predictions co-vary with the true CATE. Therefore, the p-value on the estimate of the coefficient also acts as an omnibus test for the presence of heterogeneity. If the coefficient is significantly greater than zero, then we can reject the null hypothesis of no heterogeneity.

However, coefficients smaller than 0 are not meaningful and should not be interpreted (source).

GATES

We are interested in the estimation of the GATES for two main reasons:

• They quantify the extent to which the effects differ across groups;

• They allow to assess whether we detect systematic heterogeneity or just estimation noise by testing a set of hypotheses.

If the slices are built on the true CATE, we would observe the following monotonicity:\[ \gamma_1^* \leq \ldots \leq \gamma_K^* \]where the \(^*\) indicates that the parameter builds on the true CATE. Therefore, we expect to see the same monotonicity if \(\hat{\tau}(X)\) provides a good approximation of \(\tau(X)\).

summary(evaluation, target = "GATES")

plot(evaluation, target = "GATES")

When there isn’t much detectable heterogeneity, the plot above can end up being non-monotonic. This can mean that the number of observations is too small for us to be able to detect subgroups with relevant differences in treatment effect.

In our case, the favorable monotonicity is in question. For most estimators, we can reject the null hypothesis \(H_0: GATES_5 = GATES_1\) at the 10% significance level, but not at the 5% level.

TOC curve

We are interested in the presence of heterogeneity and how much benefit there is to prioritizing treatment based on this heterogeneity. By “benefit” one means the expected increase in outcomes from giving treatment to a fraction of the population with the \(\hat{\tau}(x)\) as opposed to giving treatment to a randomly selected fraction of the same size.

The TOC is defined as the Total Outcome Contrast, expressed as:

\[ TOC(u; \hat{\tau}) = E[Y(1) - Y(0) | F(\hat{\tau}(X)) \geq 1 - u] - E[Y(1) - Y(0)] \]

This measure contrasts the average effect in the subgroup with the top 100u% highest estimated Conditional Average Treatment Effects (CATEs) with the Average Treatment Effect (ATE).

If the estimated CATEs effectively capture systematic heterogeneity, we anticipate:

• the TOC to be positive throughout, indicating that the average effect in the subgroup with high estimated CATEs consistently exceeds the overall average treatment effect.

• the TOC to be largest for low values of u, signifying that the contrast in average effects is most pronounced when focusing on the subgroup with the highest estimated CATEs.

plot(evaluation, target = "TOC")

If there is scarcely any heterogeneity in \(\tau(X_i)\), this area will be vanishingly small. This is practically what is observed. There are barely any benefits to stratifying treatment.

A null effect can be attributed to

1. Low power, as perhaps the sample size is not large enough to detect HTEs,
2. That the HTE estimator does not detect them, or
3. The heterogeneity in the treatment effects along observable predictor variables are negligible.

I would like evaluCATE to include an option to restrict the Y-axis. In this case, the initial instability makes the analysis challenging.

Turning a blind eye to the TOC curve, the results seem to be modestly significant and even promising. Heterogeneity is uncovered, ladies and gentlemen. This assignment was truly a pleasure.

Causal Inference meme knowledge wasn’t in vain

Was any heterogeneity uncovered, though? Or is this a case of seed hacking and pure luck?

8. Results variability

To explore variability, I perform 100 “Monte Carlo” simulations—i.e., I run the entire assignment 100 times using different seeds. The parameter of interest is \(\beta_2\) in BLP. I count how many times its confidence intervals cover zero, employing both Weighted Residual and Horvitz-Thompson BLPs (with de-noising).

Additionally, I conduct the same operation for three different sample splitting rules: 50/50, 70/30, and 30/70.

Each column takes around 40 minutes to run; therefore, the code is commented.

# coverage = monte_carlo(X,W,Y,iterations = 100,train_share = 0.5)
# coverage

Table 2: Third Treatment Arm - Multiple Estimations with Varying Seeds and Train/Test Splitting

Method	Train-test: 50/50	Train-test: 70/30	Train-test: 30/70
T-learner Weight	0.11	0.15	0.13
T-learner HT	0.11	0.15	0.17
X-learner Weight	0.08	0.19	0.20
X-learner HT	0.09	0.20	0.19
Causal Forest Weight	0.13	0.17	0.13
Causal forest HT	0.14	0.18	0.10
DR-learner Weight	0.03	0.10	0.13
DR-learner HT	0.03	0.11	0.15
R-learner Weight	0.03	0.07	0.13
R-learner HT	0.04	0.10	0.14

We are able to detect heterogeneous effects only sporadically. Either the effects are homogeneous or CATE estimates are not reliable.

9. Moving the goalpost

To reiterate: the last section has shown that there is either no heterogeneity or that we are unable to detect it. Detection issues arise from either a weak signal or insufficient sample size (as we know, Machine Learning methods are data-hungry).

In the following section, to “boost power,” I would argue that the treatments are quite similar. One can combine them into “no treatment” and “some treatment,” resulting in an increased sample size: 10,559 treated and 3,354 in the control group.

The parts for the combined treatment are commented throughout the code to keep the assignment concise (and to manage the running time). Simulated results are provided in the tables.

Table 3: ATE Estimated by Naive Mean Comparison

Treatment	Estimated ATE	Std. Error	t value	P-value
Treatment 1	-0.0274	0.0386	-0.70	0.478
Treatment 2	-0.0788	0.0325	-2.42	0.0153
Treatment 3	-0.0671	0.0353	-1.89	0.0576
Treatment 4	-0.0854	0.0358	-2.38	0.0171
Treatment 5	-0.0335	0.0356	-0.94	0.3466
Treatment 6	-0.0715	0.0391	-1.82	0.0678
Combined treatment	-0.0623	0.0240	-2.58	0.0096

The idea of similar treatments is debatable, to say the least, but let’s keep it for the sake of the argument.

As an illustration, best linear predictor with combined treatment and a 70% training sample gives the following result:

Performing 100 ‘Monte Carlo’ simulations (takes the most time due to a 3-fold increased sample size) again, we observe:

# coverage = monte_carlo(X,W,Y,iterations = 100,train_share = 0.3)
# coverage

Table 4: Combined Treatment - Multiple Estimations with Varying Seeds and Train/Test Splitting

Method	Train-test: 50/50	Train-test: 70/30	Train-test: 30/70
T-learner Weight	0.05	0.03	0.04
T-learner HT	0.04	0.02	0.04
X-learner Weight	0.05	0.02	0.05
X-learner HT	0.05	0.02	0.06
Causal Forest Weight	0.05	0.02	0.07
Causal Forest HT	0.04	0.03	0.07
DR-learner Weight	0.03	0.02	0.03
DR-learner HT	0.02	0.01	0.03
R-learner Weight	0.03	0.02	0.02
R-learner HT	0.02	0.02	0.02

Here, I shall be very delicate with my interpretations. As they say in movies, “anything you say can and will be used against you in a court of law”.

Firstly, only 100 simulations were conducted. Using 95% confidence intervals, the numbers in the table closely align with the probability of a type 1 error (false positive). Personally, I have more confidence in rejecting the presence of heterogeneity.

Simultaneously, one could argue that now there is more treatment heterogeneity among different treatment arms than within a single treatment arm. Instead of adding power, the standard errors have increased, leading to a reduction in power.

In conclusion, the data have been prepared and analyzed using multiple Causal Machine Learning methods. The goal was to detect whether treatment effect heterogeneity is present.

In fact, the Average Treatment (program participation) was not significant at the 5% level, and the treatment variable also showed low predictive power for the outcome variable (unemployment spell).

Subsequently, the estimated conditional average treatment effects (CATE) were analyzed using best linear predictor, GATES, and TOC methods. Monte Carlo simulations indicate that there is no treatment effect heterogeneity (or that we were not able to detect it).

10. evaluCATE experience

To be honest, I do not have enough experience to distinguish between a good and bad R package. I didn’t experience any code issues; everything was working from the start. I would wish for more pleasing visual output (GATES looks confusing), but this is only a minor point.